% Matlab Code for replication of "Optimal Debt Maturity and Firm Investment"
% by J. Jungherr and I. Schott, September 2020
% Compute model moments

% Auxiliary variables
Y                       = PARMS.Z'.*(cap_pol2.^PARMS.PSI.*lab_pol2.^(1-PARMS.PSI)).^PARMS.ZETA;
ISS                     = PARMS.ETA*(std_pol2 + max(ltd_pol2 - debtgrids,0)).^2;
DC                      = PARMS.XI*ktilde_pol2; % default costs
I                       = PARMS.DELTA * cap_pol2;
C                       = (Y(:) - PARMS.F - ISS(:)  - DC(:) - I(:))'*probst;
L_d                     = lab_pol2(:)'*probst;
%     Equilibrium mass of firms
Mstar                   = ( PARMS.W / ( (C^PARMS.RHO) * (L_d^PARMS.THETA) ) )^(1/(PARMS.THETA+PARMS.RHO)) ;

%  Aggregate labor supply (equal to L_d * Mstar in equilibrium)
L_s                     = (PARMS.W * ( (Mstar*C)^(-PARMS.RHO)) )^(1/PARMS.THETA);
% Aggregate capital
K_d_star                = cap_pol2(:)'*probst*Mstar;
% GDP
Y_star                  = (Y(:) - PARMS.F - ISS(:)  - DC(:) )'*probst*Mstar;
% Aggregate consumption from resource constraint (equal to C * Mstar)
C_res                   = Y_star - PARMS.DELTA*K_d_star;

%% Compute distributions of key variables

% a) Leverage
var                     = lev_pol2(:);
[fsort, order]          = sort(var,'ascend');
firms                   = probst(order) / sum(probst); % Apply order to probst
idx_p25                 = find(cumsum(firms)>=0.25,1,'first');
idx_p50                 = find(cumsum(firms)>=0.5,1,'first');
idx_p75                 = find(cumsum(firms)>=0.75,1,'first');
modelmoms.leverage.mean = var'*probst / sum(probst);
modelmoms.leverage.p25  = fsort(idx_p25);
modelmoms.leverage.p50  = fsort(idx_p50);
modelmoms.leverage.p75  = fsort(idx_p75);

% b) Long-term debt share
totaldebt               = (1+PARMS.C)/(1+PARMS.R)*std_pol2 + ...
                            (PARMS.GAMMA+PARMS.C)/(PARMS.GAMMA+(PARMS.R))*ltd_pol2;
ltd_share               = 1./totaldebt(:).*(PARMS.GAMMA+PARMS.C)....
    /(PARMS.GAMMA+PARMS.R).*((1-PARMS.GAMMA)/(1+PARMS.R))^4.*ltd_pol2(:);
[fsort, order]          = sort(ltd_share,'ascend');
firms                   = probst(order) / sum(probst); % Apply order to probst
idx_p25                 = find(cumsum(firms)>=0.25,1,'first');
idx_p50                 = find(cumsum(firms)>=0.5,1,'first');
idx_p75                 = find(cumsum(firms)>=0.75,1,'first');
modelmoms.ltd_share.mean= ltd_share'*probst / sum(probst);
modelmoms.ltd_share.p25 = fsort(idx_p25);
modelmoms.ltd_share.p50 = fsort(idx_p50);
modelmoms.ltd_share.p75 = fsort(idx_p75);

% c) Long-term bond spreads
ltspr                   = ((PARMS.C+PARMS.GAMMA)./p_pol2(:) + 1 - PARMS.GAMMA).^4 - (1 + PARMS.R)^4;
stspr                   = ((1+PARMS.C)./ps_pol2(:)).^4 -  (1 + PARMS.R)^4;
var                     = ltspr(:);
[fsort, order]          = sort(var,'ascend');
firms                   = probst(order) / sum(probst); % Apply order to probst
idx_p25                 = find(cumsum(firms)>=0.25,1,'first');
idx_p50                 = find(cumsum(firms)>=0.5,1,'first');
idx_p75                 = find(cumsum(firms)>=0.75,1,'first');
modelmoms.ltd_spread.mean   = var'*probst / sum(probst);
modelmoms.ltd_spread.p25= fsort(idx_p25);
modelmoms.ltd_spread.p50= fsort(idx_p50);
modelmoms.ltd_spread.p75= fsort(idx_p75);


%% Generate firm panel

% Simulate firm panel to obtain within-firm investment rates as well as size and age bins
compute_firm_panel

% Collect moments from firm panel
modelmoms.irate_sd_med      = median(ir_stds,'omitnan'); % median of firm-level sd of net inv rate
modelmoms.i_med_net         = median(ir_mean_net,'omitnan');    % median of firm-level mean of net inv rate

% Export size bins from simulated firm panel to Stata (Figures 3,4,5,6)
panel_sizegroups
% Export age bins from simulated firm panel to Stata (Figures 3,4,5,6)
panel_agegroups

%% Collect additional model moments

% Default rate
moment.defrate      = def_pol2(:)'*probst / sum(probst);
% Total exit rate (default and exogenous exit)
moment.xrate        = dist(1,1)/sum(dist(:));
% Mass of firms
moment.mass         = sum(probst)*Mstar;
% Issuance-weighted credit spread
std_iss             = std_pol2(:); % short-term debt issuance
ltd_iss             = ltd_pol2(:) - debtgrids(:); % long-term debt issuance
iwspread            = ltd_iss./(ltd_iss + std_iss) .* ltspr + ...
    std_iss./(std_iss + ltd_iss) .* stspr;
moment.iwspread     = iwspread'*probst / sum(probst);

% Value of entry
moment.valentry     = value_pol(1);
% Average capital
moment.AvgK         = cap_pol2(:)'*probst/sum(probst);
% Aggregate capital
moment.AggK         = cap_pol2(:)'*probst*Mstar;
% GDP
moment.GDP          = Y_star;
% Consumption
moment.Cons         = C_res;

